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Abstract 

We estimate the rate of aftershocks triggered by a heterogeneous stress change, using the 
rate-and-state model of Dieterich [1994]. We show that an exponential stress distribution 
P t (t) ~ exp(— t/tq) gives an Omori law decay of aftershocks with time ~ with an exponent 
p = 1 — Aa n /ro, where A is a parameter of the rate-and-state friction law, and a n the normal 
stress. Omori exponent p thus decreases if the stress "heterogeneity" r decreases. We also invert 
the stress distribution P t {t) from the seismicity rate R(t), assuming that the stress does not 
change with time. We apply this method to a synthetic stress map, using the (modified) scale 
invariant " k 2 " slip model [Herrero and Bernard, 1994]. We generate synthetic aftershock 
catalogs from this stress change. The seismicity rate on the rupture area shows a huge increase 
at short times, even if the stress decreases on average. Aftershocks are clustered in the regions of 
low slip, but the spatial distribution is more diffuse than for a simple slip dislocation. Because 
the stress field is very heterogeneous, there are many patches of positive stress changes 
everywhere on the fault. This stochastic slip model gives a Gaussian stress distribution, but 
nevertheless produces an aftershock rate which is very close to Omori's law, with an effective 
p < 1, which increases slowly with time. We obtain a good estimation of the stress distribution 
for realistic catalogs, when we constrain the shape of the distribution. However, there are 
probably other factors which also affect the temporal decay of aftershocks with time. In 
particular, heterogeneity of Aa n can also modify the parameters p and c of Omori's law. Finally, 
we show that stress shadows are very difficult to observe in a heterogeneous stress context. 
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Introduction 

Much progress has been made in describing earth- 
quake behavior based on the predictions of rate-and- 
state friction. The rate-and-state model explains the 
1/t decay of aftershock rate as a function of the time 
t since the mainshock (Omori's law) independent of 
the mainshock magnitude, the scaling of aftershock 
duration with stressing rate, the slow diffusion of af- 
tershocks with time [Dieterich, 1994]. This success led 
several authors to provide time-dependent earthquake 
probabilities using this model [Toda et al., 1998; 2003; 
2005]. Many other physical mechanisms have been 
proposed to explain Omori law, such as sub-critical 
crack growth [Das and Scholz, 1981; Shaw, 1993], vis- 
cous relaxation [Mikumo and Miyatake, 1979], static 
fatigue [Scholz, 1968; Narteau et al., 2002], postseis- 
mic slip [Schaff et al, 1998], or pore fluid flow [Nur 
and Booker, 1972]. The rate-and-state model of Di- 
eterich [1994] is probably the best candidate, however, 
because it only relies on a rate-and-state dependent 
friction law observed in laboratory experiments. 

At the same time, a number of fundamental puzzles 
remain. One of the most striking is the abundance 
of aftershocks on the rupture surface, where indeed 
most aftershocks occur. This is in stark contrast with 
simple pictures of the rupture process, which suggest 
stress should have decreased on the rupture surface 
and there should therefore be a dearth of aftershocks 
there. A second fundamental puzzle concerns the time 
dependence of aftershocks. Here, subtle but signifi- 
cant deviations from the pure Omori law inverse time 
decay of the rate of aftershocks is seen in averages of 
aftershock rates [Helmstetter et al., 2005]. While Di- 
eterich [1994] explained this as a consequence of the 
spatial dependence of stress as it decreases away from 
the fault, or as a change of stressing rate with time, 
such mechanisms do not seem to properly explain the 
aftershocks occurring on the rupture area. Thus, both 
the spatial and temporal distribution of the majority 
of aftershocks have yet to be fully explained. Here, we 
show how an extension of the rate-and-state formu- 
lation, which takes as its foundation a heterogeneous 
stress field, can explain these observations. We then 
use this model to estimate stress heterogeneity from 
aftershock rates. 

Our work builds off of the pioneering work of Di- 
eterich [1994], who derived a relation between seis- 
micity rate and stress history, for a population of 
faults obeying rate-and-state friction. For a uniform 
positive stress step (e.g., a mainshock), the rate-and- 



state model gives an Omori law decay of the seis- 
micity rate R(t) ~ t~ p with p = 1 for intermediate 
times. At very short times, smaller than a charac- 
teristic time c, which depends on the stress change, 
the seismicity rate is constant. Dieterich [1994] also 
computed the aftershock rate for a dislocation, with 
a uniform stress decrease on the rupture area, and a 
positive stress change outside the rupture, decaying 
as r <~ l/\/f in the near field, and r ~ 1/r 3 in the far 
field for r » L. As distance from the fault increases, 
the characteristic time c (typical time between main- 
shock and aftershocks) increases. Integrating over the 
fault, the seismicity rate approximately obeys Omori 
law R(t) ~ l/£ p , with an apparent exponent p < 1. 

Dieterich et al. [2000, 2003] used the rate-and- 
state model of seismicity to invert stress history from 
seismicity rate, and apply this method to Hawaii seis- 
micity. They discretize the space, with a grid size of 
about 1 km, and assume that the stress is uniform in 
each cell. This method then gives the stress history 
in each cell. The assumption that the stress is uni- 
form at scales of a few km is reasonable for the stress 
change induced by a dyke intrusion, as in [Dieterich, 
2000, 2003], or for the coseismic stress change induced 
by a large earthquake in the far field. However, the 
coseismic stress change on the mainshock fault plane, 
where most aftershocks occur, is probably very het- 
erogeneous at all scales [Herrero and Bernard, 1994]. 

In this paper, we investigate how heterogeneity of 
the Coulomb stress change and of the normal stress 
modifies the temporal decay of aftershocks with time, 
both on the fault and off-the fault. We assume that 
the stress changes instantaneously after the main- 
shock, and we neglect the relaxation of stress on the 
fault due to aseismic slip or viscous relaxation. We 
also neglect the stress change, and seismicity rate 
change, induced by aftershocks. We then try to invert 
for the stress distribution on the fault plane from the 
aftershock rate, using the rate-and-state model, and 
assuming the main source of heterogeneity is the co- 
seismic stress change. 

Relation between stress distribution 
and seismicity rate 

Dieterich [1994] derives a differential equation which 
gives the seismicity rate R(t, r) as a function of the 
stress history r(t). His model assumes an infinite pop- 
ulation of faults which obeys rate and state friction, 
with the same properties for all faults. 
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The state variable 7 is related to the stress r by 



1 



[dt-jdr] , 



(1) 



where r is the " modified" Coulomb stress change [Di- 
eterich et al., 2000], and a n is the normal stress. The 
state variable 7 is a function of the seismicity rate 
R(t,r) 

where R r is the steady state seismicity rate at the 
reference stressing rate f r . From laboratory exper- 
iments, coefficient A generally has values between 
0.005 and 0.02, for various temperature and pressure 
conditions [Dieterich, 1994]. 

Dieterich [1994] used expression (1) to derive the 
seismicity rate R(t, r) triggered by a single stress step 
r. We assume that stress rate after the stress step is 
constant dr/dt = r r , and that the seismicity rate be- 
fore the mainshock is equal to the reference seismicity 
rate R r . Using (1), the seismicity rate following the 
stress step is 



R(t,r) = 



R r 



-t/A<t„ 



l) e-*/*« + 1 ' 



(3) 



where t a is the duration of the aftershock sequence 

Aa n 



t a 



(4) 



This relation (3) is illustrated in Figure 1 for dif- 
ferent values of the stress change. For each posi- 
tive stress value, the seismicity rate is constant for 
t <C t a e~ T / Aan , and then decreases with time for 
t a e~ T / Aan « i « i„ according to Omori law with 
an exponent p = 1. For a negative stress change, 
the seismicity rate decreases after the mainshock. In 
both cases, the seismicity rate recovers its reference 
value R = R r for t » t a . The goal of this work is 
to extract the stress distribution from the seismicity 
rate. This is a difficult problem, because, as shown in 
Figure 1, the seismicity rate does not depend on the 
stress change over a relatively large time interval. 

For a heterogeneous stress field r(r), with a distri- 
bution (probability density function) P t (t), the seis- 
micity rate integrated over space is 



R(t) 



R(t,T(r))df (5) 
R{t,T)P T {r)dT (6) 



R(t, c) P c (c) dc 



(7) 



where c = t a e T / Aa " is a characteristic time of the 
aftershock rate, such that R(t, c) <~ 1/c for t <C c and 
R(t) ~ l/t for c< t < t a . 

Equation (7) is a Fredholm integral equation of the 
first kind. It has, at most, one solution [Riele, 1985]. 
Equation (7) has a simple approximate solution in 
the case when the stress change has an exponential 
distribution 



~ e 



r/r 



(8) 



where To is a positive scaling stress parameter, which 
characterizes the width of the stress distribution. 
This corresponds to a power-law distribution of "cor- 
ner times" c 

P(c) =P t (t)^ = c -i+^W™. (9) 
dc 

We also consider an approximate expression for the 
seismicity rate (3) valid for short times i<(„ 



R(t,r) 



H, 



R r t a 



e -r/Aa n + t / ta c + f - 

Substituting (9) and (10) in (7), we get 



00 

R(t) = J 



R r t a c- 1+A(T "/ T0 



c + t 



dc 
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ll-Aa„/T 



(10) 



(11) 



Expression (11) corresponds to Omori law with an 
exponent 

P=l- — • (12) 

TO 

Because equation (6) has at most one solution, the ex- 
ponential stress distribution is the only distribution 
which produces a pure Omori law decay for (<t a , 
without any cut-off or crossover at short times. How- 
ever, other distributions, e.g., a Gaussian, produce 
aftershock rate that is very close to Omori's law, over 
a very large time range. The stress distribution for 
small or negative values is not constrained by the seis- 
micity rate at short times t <C t a , so deviations from 
an exponential for negative stresses does not produce 
deviations from Omori law at short times. 

Expression (12) shows that Omori exponent de- 
pends on stress heterogeneity. The parameter To rep- 
resents the width of the stress distribution for t > 0. 
The more heterogeneous the stress is (larger To), the 
larger p is (closer to 1). Figure 2 illustrates how the 
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rate-and-state model with a heterogeneous stress dis- 
tribution produces a power-law decay with an expo- 
nent p < 1. 

Helmstetter et al. [2005] found that, for stacked 
aftershock sequences in Southern California, Omori 
exponent is close to 0.9, for times ranging between 
a minute (but possibly even less) and one year, and 
for mainshock magnitudes between 2 and 7.5. This 
suggests that the stress distribution is close to expo- 
nential in the tail, with a characteristic stress r ~ 
10Aa n . Assuming that A = 0.01 (as measured in 
laboratory friction experiments [Dieterich, 1994]) and 
a n = 100 MPa (corresponding to the lithostatic pres- 
sure at a depth of about 5 km), this gives Aa n = 1 
MPa and t = 10 MPa, a value larger than the typical 
stress drop cxo = 3 MPa [Ide and Beroza, 2001], but of 
the same order of magnitude. However, a few studies 
tried to estimate Aa n directly from earthquake cat- 
alogs, and obtained values smaller than the ones de- 
rived from the laboratory value of A. Dieterich [1994] 
found Aa n = oo/20, from the relation between af- 
tershock duration and the recurrence time (assuming 
characteristic earthquakes). This gives Au n = 0.15 
MPa assuming a stress drop of 3 MPa. Cochran et 
al. [2004] used the rate-and-state model to model 
tidal triggering of earthquakes, and obtained a prcf- 
ered value of Aa n = 0.064 MPa, and an acceptable 
range 0.048 < Aa n < 0.11 MPa. 

The rate-and-state model with a uniform stress 
step (3) cannot explain an Omori law decay with 
p > 1. Equation (7) does not have a solution with 
P t (t) > and t <C t a in this case. Some aftershock se- 
quences however have an Omori exponent larger than 
one. The only solution in order to obtain a p-value 
larger than one in the rate-and-state model is to have 
a variation of stress with time, which may be due to 
postseismic slip or viscous relaxation, although these 
explanations involve relatively large stress changes 
with time [Dieterich, 1994]. Other explanations for 
Omori's law do allow for larger p- values [Mikumo and 
Miayatake, 1979; Shaw, 1993; Narteau et al, 2002]. 

Estimating the stress distribution from 
aftershock rate 

We have shown above that, according to the rate- 
and-state model, the Omori exponent provides some 
information on the stress heterogeneity (but only if 
p < 1). Furthermore, we can (in theory) obtain the 
complete stress distribution (in the region where we 
measure the seismicity rate) from the temporal evo- 



lution of the seismicity rate. Expression (6) indeed 
provides a method for estimating the full distribution 
P T (r), provided we observe the seismicity rate R(t) 
over a wide enough time interval. 

We first discrctize the integration over stress and 
times, using a linear sampling for stress, and a loga- 
rithmic sampling for times, using the same number N 
of points. Equation (6) is then similar to the system 
of iV linear equations 

N 

R(tj) = R(tj,n) PAn) (T i+ i - n) (13) 

i=l 

We divide both sides of equation (13) by R(tj) to 
stabilize the problem. Equation (13) thus becomes 

1 = g^^) Pr(ri)(Ti+i _ ri) = MxP (u) 

where M is a N x N matrix M(i,j) = (rj+i — 
Tj) R(tj,Ti)/R(tj) and the vector P is the stress dis- 
tribution at points t\,...,tn- 

The inversion of the stress distribution from (14) 
is an ill-posed problem, i.e., the solution is very sensi- 
tive to noise. We thus use the regularization method 
of [Riele, 1985]. We introduce an additional con- 
straint to (13), minimizing either the first derivative 
|P'(r)|, the smoothness ||P"(r)||, or the distance be- 
tween P t (t) and an initial guess Po(t). ( e -g-> a Gaus- 
sian distribution). Instead of solving directly (14), we 
minimize the quantity 

\\M P - l\\ 2 + a\\L(P)\\ 2 , (15) 

where a > is the regularization parameter, and L 
is a linear operator, e.g., L(P) = P — P , L(P) = P' 
(first derivative), or L(P) = P" (second derivative). 
We also impose that the stress distribution is positive. 
We thus search for the positive vector P that mini- 
mizes equation (15), using the non-linear least-square 
fitting program given by Lawson and Hanson [1974]. 

In practice, the estimation of P t (t) for large r 
is limited by the minimum time i m i n at which we 
can reliably estimate the seismicity rate. The largest 
stress we can resolve is of the order of r max = 
— Aa n log(t m i n /t a ). Practically, this time i m i n may 
be as low as a few seconds, if we correct from catalog 
incompleteness shortly after the mainshock [ Vidale et 
al., 2004]. For negative stress, we are limited by the 
maximum time i max after the mainshock, and by our 
assumptions that secondary aftershocks are negligi- 
ble, and that the stress does not change with time 
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(e.g., neglecting post-seismic relaxation). In order to 
resolve P t {t) for negative values, we need to know the 
seismicity rate for times larger than the aftershock du- 
ration t a (i.e., usually at least a few years). Indeed, 
the seismicity rate after a stress decrease is close to 
zero for t <C t a , so that the measure of R(t) for t <C t a 
does not provide any information on P t (t) for t < 0. 

Application of the method to a 
stochastic slip model 

Stochastic k 2 slip model 

We have tested the rate-and-state model on a re- 
alistic synthetic slip pattern. Herrero and Bernard 
[1994] proposed a kinematic, self-similar model of 
earthquakes. They assumed that the slip distribu- 
tion at small scales, compared to the rupture length 
L, does not depend on L. This led to a slip power- 
spectrum for high wave-number equal to 

u(k) = C—^r forfc>l/L, (16) 
ji k z 

where a is the stress drop (typically 3 MPa), [i is 
the rigidity (typically 3300 MPa in the lower crust), 
and C is a shape factor close to 1. For wavelengths 
larger than the rupture length L, the power spectrum 
is constant 

u{k) = C— L 3 forfc<l/i. (17) 
I- 1 

This model (16) reproduces the l// 2 power-spectrum 
of seismograms for large frequencies [Herrero and 
Bernard, 1994]. 

Shear stress change and seismicity rate on 
the fault 

We have used the k 2 model to generate a synthetic 
slip pattern, and compute the shear stress change on 
the fault from the slip [Andrews, 1980; Ripperger and 
Mai, 2004]. Note that the seismicity rate given by 
(3) depends on the Coulomb stress change, which is 
equal to the shear stress change on the fault because 
the normal stress change on a planar fault is zero. 
If we analyze off-fault aftershocks or complex rupture 
geometries, we would have to consider changes in nor- 
mal stress as well. 

We have modified the k 2 model in order to have 
a finite standard deviation of the stress distribution. 
The k 2 model (16) produces a shear stress change 
with a power spectrum r(fc) <~ k^ 1 for large k, be- 
cause the stress is approximately the derivative of the 



slip. As a consequence, the shear stress change for 
the k 2 model is extremely heterogeneous, with an in- 
finite standard deviation. The exponent n = 2 in the 
k 2 model (16) is thus a minimum physical value for 
the slip power-spectrum [Herrero and Bernard, 1994]. 
Using u(k) ~ k~ 2 produces a shear stress change with 
a standard deviation which diverges logarithmically 
as the maximum wavenumber increases. Thus Omori 
p- value for this slip model tends to 1 as the grid reso- 
lution increases. We have thus replaced the exponent 
n = 2 in (16) by n = 2.3, and smoothed the crossover 
at k — l/L, using 

u (k)=C^ nT L \ , ■ (18) 
w /x (kL + 1)" v ' 

We have computed the stress change on the fault 
from this synthetic slip model, for a fault of 50 x 50 
km, with a resolution dx — 0.1 km, and a stress drop 
(To = 3 MPa (i.e., the average stress change on the 
fault is —3 MPa). The maps of the slip and stress on 
the fault arc shown in Figure 3. The stress field has 
large variations, from about -90 to 90 MPa, due to 
slip variability. We did not constrain the slip to be 
positive. This could be done by changing the phase 
of the lowest mode, and tapering the slip close to the 
edges, so that the maximum slip is at the center [Her- 
rero and Bernard, 1994]. Doing so introduces small 
deviations of the stress distribution from a Gaussian 
distribution for r m 0, but does not introduce signifi- 
cant changes on the seismicity rate . 

We have then estimated the seismicity rate on the 
fault predicted by the rate-and-state model, by inte- 
grating numerically (5) using the observed stress map, 
and Aa n = 1 MPa. While the stress on average de- 
creases on the fault, the seismicity rate shows a huge 
increase after the mainshock (by a factor 10 10 , but, 
of course, the seismicity rate at short times, smaller 
than the duration of the earthquake, has no physi- 
cal sense) (see Figure 4). It then decays with time 
approximately according to Omori law, with an ap- 
parent exponent p = 0.93. At large times t w t a , the 
seismicity rate decreases below its reference rate due 
to the negative stress values. 

Marsan [2006] reached similar conclusions, using 
the same model: the main effect of stress heterogene- 
ity on the fault is to produce a short term increase 
of the seismicity rate, and to delay the seismic quies- 
cence on the fault by months to years. 
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Synthetic aftershock catalog 

We have generated synthetic earthquake catalogs 
according to the rate-and-state model, using the (mod- 
ified) k 2 model (18) to generate the stress change. 
We have simulated aftershock sequences triggered by 
this heterogeneous stress change, without including 
earthquakes interaction (i.e., without coseismic stress 
changes induced by aftershocks), using the method of 
Dieterich et al. [2003]. We assume a non-stationary 
Poisson process with an average seismicity rate R(t, r) 
given by (3). We generate aftershock in each cell inde- 
pendently of the other cells, assuming that the stress 
is uniform in each cell. We do not need to gener- 
ate event magnitudes, because we do not include sec- 
ondary aftershocks in our simulation. We consider 
that each aftershock does not modify the stress field 
or the seismicity rate. 

In each cell, we generate events one after the other. 
If the last event in the cell occurred at a time tj after 
the mainshock , the probability that the next earth- 
quake will occur at a time smaller than ti + dt is given 

by 



F(dt, U) — exp 



ti+dt 



R(t',r) dt' 



(19) 



The function F{dt ) t i ) increases from to 1 as dt 
ranges from to oo. To determine the time t i+i = 
ti+dt of the next event, we generate a random number 
z between and 1, and we solve for F(dt, ti) — z. We 
have generated 6 synthetic catalogs from the stress 
field shown in Figure 3, using Aa n = 1 MPa or 
Aa n = 0.1 MPa. We used different values values of 
the reference rate R r , and of time interval t m i n — t max 
(see Table 1), in order to test how the inversion 
method depends on the quality of the catalog. 

Inversion of stress history from seismicity 
rate 

We have first applied the method of Dieterich et 
al. [2000, 2003] on this synthetic stress field shown 
in Figure 3b. Dieterich et al. [2000, 2003] estimate 
the stress history r{t) at any point on a grid, assum- 
ing that the stress change is homogeneous in each 
cell, but may change with time. The stress history is 
obtained from the seismicity rate by solving equation 
(1). We wanted to apply this method on this synthetic 
stress model to test how stress heterogeneity affects 
the inverted stress change. The results are shown in 
Figure 5. The inverted stress change at short times 
is close to the maximum stress change « 100 MPa, 



and then decreases down to a value close to the av- 
erage stress change ~ —3 MPa at large times t > t a . 
Dividing the fault into smaller size cells would not im- 
prove the results very much. Because this slip model 
is self-similar, there are almost everywhere some parts 
of the fault where the stress (and thus the seismicity 
rate) increases. This shows that a small-scale stress 
heterogeneity, without any time dependence, is inter- 
preted by this method as a variation of stress with 
time. Also, it shows that a stress decrease cannot be 
resolved if it is mixed with a stress increase, unless 
looking at very long times. This may explain why 
stress shadows are so difficult to observe [Felzer et 
al, 2005]. 

Inversion of stress distribution from 
seismicity rate 

This test shows that variability with time is hard 
to distinguish from small-scale heterogeneity in space 
based on the temporal evolution of the seismicity rate. 

In order to characterize the coseismic stress change 
on the fault plane, we thus need to neglect one effect 
(small-scale heterogeneity) or the other (time vari- 
ation). Our method estimates the stress distribu- 
tion on the fault from the seismicity rate, assuming 
that stress does not change with time. In theory (if 
we had an infinite time interval, a huge number of 
aftershocks, no foreshocks or secondary aftershocks, 
and if we knew the parameters R r , t ai and Acr„), 
this method provides the distribution of stress on the 
fault. If the fault is divided into smaller cells, this 
method gives a map of the average stress change in 
each cell, as well as its variability. 

For each synthetic catalog, we have measured the 
seismicity rate on the fault by smoothing aftershock 
times. We used a kernel method to estimate R(t) 
from aftershocks time ti, with i = 1 to N, with a 
log-normal filter 

JV 



^) = E7T^ex P| - 



(logioffl-logiofe)) 2 
2K 2 

(20) 

with a kernel width h = 0.08. 

We then used the inversion method described pre- 
viously to estimate the stress distribution P T (r) from 
the seismicity rate. We used the regularization condi- 
tion L(P) — P' in (15), i.e., minimizing the derivative 
of P t (t), using a = 10 4 (decreasing a produces huge 
fluctuations of P t (t)). 

We have also estimated the Gaussian stress distri- 
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bution that best fits the observed seismicity rate. We 
evaluate the mean — <to, and the standard deviation 
t* of the Gaussian function, as well as the aftershock 
duration t a , using a maximum likelihood approach. 
We maximize the log-likelihood function defined by 



£ma: 

J2 log aft)- / 

— 1 AT J 



R(t)dt , 



(21) 



i=l,N 



where the seismicity rate is given by 



e -(T+a„) 2 /2r* 2 

R(t) = I R(t, t) = dr . (22) 



The log-likelihood function is maximized when the 
rate estimate R, weighted logarithmically, is large 
when events occur at times U, and when the total 
rate estimate integrated over time is low. 

Table 1 gives the parameters of each simulation, 
and the results of the inversion. Figures 7 and 9 
show the real stress distribution (evaluated from the 
stress map shown in Figure 3b) and the inverted one, 
for each synthetic aftershock catalog. We test both 
inversion methods, either solving (15) for P t (t) for 
80 < t < 80 MPa, or assuming a Gaussian stress 
distribution. 

Figure 4 compares the theoretical seismicity rate 
given by (5) using the observed stress field, with the 
seismicity rate estimated from the seismicity catalog 
using (20), and with the reconstructed seismicity rate 
estimated using (6) from the inverted stress distribu- 
tion. For this synthetic catalog, the seismicity rate is 
almost indistinguishable from an Omori law with an 
exponent p = 0.93 for t/t a < 10. 

In the first 2 catalogs in Table 1, with more than 
several thousands events, we obtain a very good es- 
timation (error less than 6%) on all parameters t* , 
(To and t a . If the number of events decreases to 392 
events, without changing the time interval, we still 
obtain a rather good estimation of t a and t* , but the 
error on the stress drop increases (see model 3 in Ta- 
ble 1). For a shorter catalog (#4 in Table 1), with 
292 events and only 4 ranges of magnitude in time, 
the stress drop is not constrained, unless we fix the 
aftershock duration to its true value. Alternatively, 
we can fix the stress drop and obtain a rather good 
estimation of r* and t a . This shows that the main ef- 
fect in recovering cr , T * and t a is the catalogue time 
interval, which needs to extend over a reasonable frac- 
tion of t a - This is because very different values of t a 



and do can produce very similar seismicity rate R(t) 
for t < t a /100, as can be shown in Figure 8. If we 
decrease Aa n , keeping r* fixed, the Omori exponent 
becomes closer to 1, and the error on all parameters 
increases (see models #5 and #6). 

When inverting for the complete distribution P t (t), 
the results are pretty good for the first simulation, 
with an unrealistic large time interval and number of 
events. There are however deviations in the tails, for 
r > 30 MPa, which correspond to very short corner 
times c — t a exp(— r/Aa n ) = 10~ 13 , much smaller 
than the minimum time t min /t a — 10 -10 used for the 
inversion of P t (t). For catalogs #2-4 in Table 1, the 
distribution of P t (t) is not constrained for t < 0, 
and for r ^> 1, because of the limited time interval. 
The results are very poor for both simulations #5 
and #6 in Table 1, with Aa n = 0.1 MPa and Omori 
exponent p = 0.993. In this case, we have almost 
no resolution on P t (t) for t < 0. This method only 
provides a rough estimate of the width of the distri- 
bution for r > 0. Thus, in practice, unless one has 
a very long catalogue in time, and significant devi- 
ations from Omori law, little can be said about the 
stress shadow regions. 

Off-fault aftershocks 

We can make simple estimates of the stress change 
and seismicity rate off of the fault plane. For mode 
III rupture, static elasticity reduces to a Laplacian 
Am = 0. For a Laplacian, a Fourier mode with 
wavenumber k along an infinite fault decays expo- 
nentially into the bulk proportional to k times the 
distance y to the fault. With these basis functions, 
we can easily extrapolate off of the fault, although 
since it neglects rupture end effects, it is valid only 
for distances less than the rupture length L and in 
areas along-side the mainshock rupture area, and not 
extending into the lobes of increasing stress beyond 
the finite rupture length. Thus, we are looking at 
regions which would be in the stress "shadow" of a 
simple rupture. Within this region, at a distance 
y < L from the fault, the power-spectrum of the dis- 
placement for the modified slip model (18) becomes 
u(k,y) ~ (kL + l)~ n exp(—ky). The power-spectrum 
of the stress change is given by 



r(k,y) ~ ku{k,y) 



-ky 



(fcL+1) 



n-1 ' 



(23) 



This shows how the stress heterogeneity decays very 
rapidly with distance from the rupture surface. 
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Figure 10a shows the seismicity rate for different 
values of the distance from the fault y/L, using the 
slip model shown in Figure 3a. We computed the 
stress at a distance y from the fault using r(k, y) = 
r(fc, 0) exp(— k y), i.e., multiplying the stress map 
shown in Figure 3b by exp(— ky) in the Fourier do- 
main. The stress distribution is reasonably close to 
a Gaussian distribution at all distances. Therefore, 
we have used the best-fitting Gaussian distribution 
in order to compute the seismicity rate shown in Fig- 
ure 10a. The standard deviation of the stress distribu- 
tion decreases very fast with the distance to the fault, 
which produces a strong drop of the seismicity rate off 
of the fault. The average stress decreases much slower 
with y. Figure 10b shows the falloff with distance of 
these quantities. For y/L > 0.1, the stress field is 
much more homogeneous and mostly negative (the 
standard deviation is smaller than the absolute mean 
stress). Therefore, the seismicity rate for y/L > 0.1 
is smaller than the reference rate at all times t < t a . 
Marsan [2006] also used the rate-and-state model to 
investigate how stress change heterogeneity modifies 
the rate of off-fault aftershocks. He used a slightly 
different slip model, and assumed the spatial depen- 
dence of the stress variability decayed with the same 
form as the stress, as distance cubed. With this as- 
sumption he found, not surprisingly, larger distances 
of triggering. 

In practice, it is difficult to analyze the rate of 
off-fault aftershocks, because the aftershock rate and 
the reference seismicity rate decrease with the dis- 
tance from the fault, and because of location er- 
rors. Also, secondary aftershocks triggered by off- 
fault events will perturb the stress field and seismicity 
rate with additional stress heterogeneity. Our seis- 
micity rate estimates here presume focal mechanisms 
of aftershocks similar to the mainshock focal mech- 
anism; other focal mechanisms could have different 
rates, but optimally oriented plane estimates may not 
be the best approach [McCloskey et at, 2003]. In any 
case, we do see very rapid falloff of the seismicity with 
distance from the fault, a point which deserves fur- 
ther observational exploration. Note that Figure 10a 
shows the seismicity rate normalized by the reference 
rate R r . If R r decreases with the distance to the fault, 
the decrease of the aftershock rate with y will be even 
faster than shown in Figure 10a. 



Gaussian versus exponential stress 
distribution 

While the pure Omori law with p < 1 occurs for the 
exponential distribution of stress changes, we find nu- 
merically that a Gaussian stress distribution (which 
the k 2 model and many other models give) , also gives 
realistic looking p-values over wide ranges of time 
scales. Some insight into why this is the case can 
be gained by noting that for a sufficiently wide range 
of values, a Gaussian is a good enough approximation 
of an exponential. Taking the ratio of a Gaussian to 
an exponential 



exp 



cxp 



(t + ap) 1 
2t* 2 

1 (t + (t 



/exp — 

T 



T 
70 



CTo 
TO 



2t 2 



(24) 



For t/t* = t*/t - (to/t* ± 1 this is within a fac- 
tor exp(l) of being constant. Thus, over an e-folding 
range of t* /tq we have something well approximated 
by an exponential. 

We can use this result to obtain an approximate 
analytical expression for the effective Omori expo- 
nent produced by a Gaussian stress distribution. Ex- 
pression (24) shows that the exponential distribution 
closer to the Gaussian one for a stress r has a char- 
acteristic parameter t — T* 2 /(a + t) . As Figure 2 
illustrates, the more important contribution to the af- 
tershock rate at time t <C t a is due to stress values of 
the order of t c = — Aa n log(t/t a ). If the stress change 
obeys a Gaussian distribution, stresses larger than t c 
are less frequent than for r = t c , therefore they have 
a smaller contribution to the seismicity rate at time 
t. Smaller stress values t < r c are more frequent, but 
the seismicity rate at time t is negligible compared to 
larger stress values. We thus obtain the following rela- 
tion between the parameters (T and r* of a Gaussian 
distribution, and the parameter To of the exponential 
distribution which better explains the aftershock rate 
at a given time t 



To = 



(To - Aa n \og(t/t a ) 



(25) 



Using expression (12), we obtain the following rela- 
tion between the effective Omori exponent at time t 
and the parameters Co and t* of a Gaussian stress 
distribution 



1 - 



Aa n a - A 2 a 2 n log(t/t a ) 



(26) 
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showing the slow increase of p with time. Figure 11 
compares this approximate solution (26) with the 
variation of p with time computed by integrating nu- 
merically (6), using Aa n = 1 MPa, for a Gaussian 
stress distribution with (T = 3 MPa and r* = 10 
MPa. The approximate solution (26) for Omori ex- 
ponent is quite good for short times t -C t a , but the 
difference with the exact solution increases as time 
approaches t a - This expression (26) also shows us the 
inherent tradeoff between the mean stress change —do 
and the variance of the stress change r* in affecting 
the p- value. 

Discussion 

We have considered above only heterogeneity of the 
Coulomb stress change. However, there are other im- 
portant factors that affect the temporal evolution of 
the seismicity rate, such as heterogeneity of the fric- 
tion law parameter A, normal stress, and stressing 
rate, multiple interactions between aftershocks, fore- 
shocks, and postseismic relaxation. 

Heterogeneity of the friction parameter A, 
normal stress, stressing rate, and reference 
seismicity rate. 

We have shown that Coulomb stress change het- 
eroegeneity modifies the temporal evolution of the 
seismicity rate, compared to a uniform stress change. 
But other kinds of heterogeneity may also impact the 
aftershock decay with time, in particular the normal 
stress. Normal stress heterogeneity enters the prob- 
lem in two ways, through the "modified" Coulomb 
stress change r, and through the Aa n term in the 
denominator. 

Slip on a rough fault will produce coseismic changes 
of the normal stress [Dieterich, 2005]. For coseismic 
changes of the normal stress which are small com- 
pared to the normal stress, we can assume that Aa n 
docs not change with time, and account for coseismic 
changes of <r n only in the coseismic Coulomb stress 
change r. For larger coseismic changes of normal 
stress, we have to use a more complex form for the 
relation (1) between stress history and seismicity rate 
[Dieterich, 1994], and equation (3) is no more valid. 

In addition to coseismic stress changes of cr„, there 
are also spatial fluctuations of Aa n . For instance, we 
expect both A and a n to change with depth. With 
a wide variety of materials making up fault zones 
and the presence of fluids, there is probably no lower 
bound on Aa n . The first effect of introducing hetero- 



geneity of Aa n is to increase the fluctuations of the 
normalized stress r/Aa n , i.e., the standard deviation 
of r/Aa n is larger than T*/(Aa n ). Therefore, Omori 
exponent increases if Aa n is more heterogeneous. Ne- 
glecting the fluctuations of Aa n will thus overestimate 

T*. 

The second effect is to introduce fluctuations of the 
aftershock duration t a , which scales with Aa n . Fluc- 
tuation of Aa n thus decrease the apparent aftershock 
duration. The value of t ai inverted assuming r/A(j n 
is uniform, is smaller than its average value t a . Also, 
Aa n heterogeneity widens the duration of the quies- 
cence (time period when R(t) < R r ). 

We have illustrated the effect of normal stress het- 
erogeneity in Figure 12, which compares the seismic- 
ity rate with and without fluctuations of Aa n . Fluctu- 
ations of coseismic Coulomb stress changes are mod- 
eled by a Gaussian distribution of mean — a = —3 
MPa and standard deviation r* = 5 MPa. For Aa n , 
we use a lognormal distribution of mean (Aa n ) = 1 
MPa and standard deviation std(Aa n ) = 7.3 MPa. 
The main effects of Aa n heterogeneity is to increase 
the apparent Omori exponent (measured for t < 0.01 
yr) from 0.44 to 0.91, and to decrease the apparent 
aftershock duration, (defined as the time when the 
aftershock rate decreases below its reference value) 
from 0.2 to 0.05 yr. 

Inverting for the Coulomb stress distribution from 
R(t), assuming that Aa n = 1 MPa everywhere, gives 
t* = 15.7 MPa, cr ,9 = 9.9 MPa and t a , s = 0.15 yr, 
instead of the true value r* = 5 MPa, a^, g = 3 MPa 
and (t a ) = 1 yr. The errors in the inverted parame- 
ters Tg , CTo, g , t a ^ g are negligeable when the Coulomb 
stress change r is more heterogeneous than Aa n , i.e., 
if r* 3> std(Aa n ). The fit of the aftershock rate with 
Aa n assumed constant gives a reasonably good fit to 
the seismicity rate computed including Aa n hetero- 
geneity. The misfit will probably be within the noise 
level for real data set. This shows that, with the time 
dependence of the seismicity alone being the source 
of information, we cannot distinguish between hetero- 
geneity of t or Aa n . Finding other effects which might 
be able to separate out these contributions of shear 
stress heterogeneity and normal stress heterogeneity 
remains an area worthy of further inquiry. 

The fact that Aa n heterogeneity increases the 
Omori exponent may explain why very low p-values 
are seldom observed, even outside the mainshock rup- 
ture area, where Coulomb stress change is relatively 
uniform (sec Figure (10b). This also explains why 
the crossover time c for off-fault aftershocks is usually 
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very short, as for on-fault aftershocks. It also makes 
stress shadows even more difficult to detect. Even 
in the regions where stress change is negative and not 
too heterogeneous (r* < cr ), fluctuations of Aa n pro- 
duce an increase of the seismicity rate at short times, 
while a uniform value of Aa n gives a quiescence at all 
times. 

Another parameter that affects the seismicity rate 
is the reference stressing rate, which modifies the af- 
tershock duration t a . Heterogeneity of the stressing 
rate will thus also yield an error in the inverted val- 
ues of t* and do- In contrast, the seismicity rate does 
not depend on the spatial fluctuations of the refer- 
ence seismicity rate, but depends only on the average 
value of R r . In practice, R r is measured from the 
average seismicity rate over a long time period before 
the mainshock. The uncertainty on R r is generally of 
a factor of about 2. This could induce large relative 
errors on the stress drop estimate ao, g , but does not 
affect too much the inverted values of r* and t a . 

Foreshocks 

An assumption of our model is that the seismicity 
rate before the mainshock is equal to the reference 
seismicity rate. But most mainshocks are preceded by 
foreshocks, so that the seismicity rate Rq before the 
mainshock is usually larger than the reference rate 
R r . Using the results of Dieterich [1994], we can take 
into account this effect by replacing the term er T l Aa ^ 
in (3) by 

l cxp (^)^ cxp (-i:- log (i;))- 

(27) 

The effect of increasing Rq is thus equivalent to shift- 
ing the stress distribution toward larger values, by 
the amount Aa n \og(R /R r )- Not correcting for this 
effect will thus over-estimate the stress change. 

Secondary aftershocks 

We have neglected in this study the role of after- 
shocks in changing the seismicity rate and redistribut- 
ing the stress. We know that most aftershocks may 
be secondary aftershocks, triggered by previous af- 
tershocks [Felzer et al., 2003; Helmstetter and Sor- 
nette, 2003]. Ziv and Rubin [2003] studied a quasi- 
static fault model that is governed by rate- and state- 
dependent friction. They have shown that, if the 
mainshock is modeled as a uniform stress increase, 
the effect main of secondary aftershocks in the rate- 
and-state model is to renormalize the seismicity rate 



without changing its time dependence (i.e., without 
changing Omori p value). If the stress change induced 
by the mainshock is non-uniform, multiple interac- 
tions between earthquakes modify the spatial distri- 
bution of aftershocks [Ziv, 2003]. 

Marsan [2006] also performed numerical simula- 
tions to model the effect of multiple interactions. He 
modeled the stress change induced by each aftershock 
by a Gaussian white noise of zero mean, i.e., assuming 
all aftershocks have the same size, and neglecting spa- 
tial correlation of the stress field. He concluded that 
the main effect of multiple interactions is to increase 
the reference rate, but also to decrease the ratio of 
the aftershock and background rates. The existence 
of multiple interactions also decreases the apparent 
aftershock duration, but does not change the Omori 
exponent. 

Therefore, secondary aftershocks should not change 
the value of the width of the stress distribution in- 
verted from the aftershock decay on the mainshock 
fault, which is controlled by Omori exponent. But 
multiple interactions may bias the value of the av- 
erage stress change and aftershock duration. Devel- 
oping more realistic models for multiple interactions 
remains an area worthy of further inquiry, but beyond 
the goals of this paper. 

Conclusion 

We have shown how a new extension of the rate- 
and-state friction formulation for seismicity rates, 
which takes as its foundation a heterogeneous stress 
field, can explain the most prevalent and puzzling of 
aftershocks, those which occur on the mainshock rup- 
ture area, where the stress decreases on average af- 
ter the mainshock. With this point of view, subtle 
but significant deviations from the pure inverse time 
omori exponent are mapped onto measures of stress 
change heterogeneity on the fault. This contrasts with 
the established methodology of Dieterich et al. [2000, 
2003] , in which these deviations are mapped onto time 
dependent stress changes following the mainshock. 

Taking the point of view that stress heterogeneity 
can be quite large at the local scale on the fault sur- 
face which ruptured, we have gained insights into a 
number of topics of relevance to stress heterogeneity's 
and earthquake behavior. Regarding stress shadows, 
we have seen how they are very difficult to detect in a 
heterogeneous stress context, relying on subtle details 
in the seismicity rates at times of order t a , subtleties 
which would become even more difficult to detect if 



11 



t a were nonuniform. 

Regarding stress heterogeneity relative to mean 
stress changes, we have found in our efforts to model 
seismicity changes with scale invariant slip distribu- 
tions that typical stress changes are actually larger 
than mean stress drops on faults, so that a picture of 
a very rough stress distribution on a fault which has 
broken is a much better picture than standard crack- 
like models, which tend to concentrate aftershocks at 
the edges of ruptures. This provides important con- 
straints on physical models of earthquakes. Finally, 
we have shown that modest catalogue lengths allow 
an accurate inversion for some stress heterogeneity 
parameters, if the only source of heterogeneity is the 
Coulomb stress change. 

However, there are probably other important fac- 
tors that affect the temporal evolution of the seis- 
micity rate, such as heterogeneity of the friction law 
parameter A, effective normal stress, and stressing 
rate, multiple interactions between aftershocks, and 
postseismic relaxation. In particular, heterogeneity of 
Aa n may explain why Omori exponent and character- 
istic time c does not seem to depend on stress change 
amplitude [Felzer, 2005]. We have shown that earth- 
quake triggering is not only controlled by the average 
values of the Coulomb stress change, or of the effec- 
tive normal stress, but rather by their heterogeneity. 
Particularly, short time aftershock rate is mainly con- 
trolled by the maximum stress change in this region, 
rather than by its average value. Estimation of coseis- 
mic slip from seismograms or geodesy is not accurate 
enough to estimate small scale variations of the stress 
change on the fault plane. Therefore, we need to cou- 
ple large scale deterministic slip models with small 
scale stochastic slip models, in order to reproduce the 
spatio-temporal distribution of triggered earthquakes. 
This also shows the difficulty of inverting the stress 
field from the spatio-temporal variations of the seis- 
micity rate. Real data is limited in number of events, 
catalog duration, and location accuracy; and we have 
only rough estimates of the large-scale average value 
of the friction parameters, normal stress and stress- 
ing rate. With the limited information given by the 
seismicity rate, it is hard to characterize the different 
factors that control earthquake triggering, especially 
on or close to the fault, where stress and material 
properties are very heterogeneous. 
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Table 1. Parameters of the synthetic aftershock catalogs: number N of events, time interval [t m i n t max ], 
value of Aa n used for the simulations, Omori exponent p (measured by maximum likelihood from the 
simulated catalog for t/t a < 0.01), average stress change — Co, and standard deviation r* (in MPa), 
and results of the inversion: oo.g, T g anc * ta,g, estimated assuming a Gaussian stress distribution P T (r). 
Stress values are in MPa. 
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a This catalog is a subset of catalog #2, obtained by increasing the minimum magnitude by one unit. 
^The parameter was fixed to its real value in the inversion 




Figure 1. Seismicity rate R(t, r) (normalized by the reference seismicity rate R r ) as a function of time (normalized 
by the aftershock duration t a ), given by the rate-and-state model with a uniform stress step (3), for different values 
of the stress change, ranging from r/Aa n = —30 (bottom) to r/Aa n — 30 (top). 
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Figure 2. The thin colored lines are the seismicity rate R(t, r) for a uniform stress change r, ranging from 
t = (blue flat curve) to r = 50 MPa (red curve), weighted by the probability -Pt(t), using Aa n = 1 MPa. The 
stress distribution is given by -FV(t) ~ exp(— r/5) with r > 0. The solid black line is the total seismicity rate 
R(t) = J °° R(t,r) P t {t) dr. The superposition of curves R(t,r) with a power-law distribution of crossover times 
c = t a cxp(— t / Aa n ) gives rise to a power law decay of R(t) with an exponent p w 0.8. The dashed lines are Omori 
laws with p = 1 (bottom) and p = 0.8 (top). 




10 


20 30 


40 




x (km) 






slip (m) 









2 4 6 


8 10 





5 


(b) 






10 




15 




20 




25 




N 






30 




35 




40 




45 





¥ * . ■ m 
' \ 




■ *. ■ ' ■ •. , • i» . *. ■ . 


L 1 




*■ 




i ._»• 


.*> . 


10 


20 30 40 
x (km) 

t (MPa) 





-50 



50 



Figure 3. (a) Stochastic slip model, with a power-spectrum u(k) = l/(kL+ 1) (where L = 50 km is the rupture 



length), a stress drop <7o 
direction). 



3 MPa, and a cell size dx = 0.1 km. (b) Shear stress change (parallel to the slip 
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Figure 4. Seismicity rate given by the rate-and-state model [Dieterich, 1994], for the stress change shown in 
Figure 3, assuming Aa n = 1 MPa, and without earthquake interactions. The solid red line is the seismicity rate 
estimated from the simulated earthquake catalog (see model #1 in Table 1). The dashed black line is a fit by 
Omori's law for t < i a /100, with exponent p — 0.93. The crosses show the fit with the rate-and-state model 
assuming a Gaussian P t (t). 
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Figure 5. Shear stress change, as a function of the time after the mainshock, estimated from the seismicity rate 
using Dieterich et al. [2000, 2003] method. We solved equation (1) for the stress history, assuming that the stress 
is uniform in space but changes with time. 
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Figure 6. Seismicity map, for a synthetic catalog generated using the rate-and-state model [Dieterich, 1994], for 
the stress change shown in Figure 3, assuming Aa n = 1 MPa, and without earthquake interactions. Only events 
with t < t a are shown. 
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Figure 7. Stress distribution estimated directly from the stress map shown in Figure 3b (circles), and inverted 
from the seismicity rate shown in Figure 4. The solid black line is the solution of equation (15)). The best fitting 
Gaussian distribution is shown as a blue dashed line. 
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Figure 8. Seismicity rate given by the rate-and-state model [Dieterich, 1994], for the stress change shown in 
Figure 3, assuming Aa n = 1 MPa, and without earthquake interactions. The solid red line is the seismicity rate 
estimated from the simulated earthquake catalog (see catalog #4 in Table 1). The crosses show the fit with the 
rate-and-state model assuming a Gaussian P T (r), and inverting for t a , <tq, and r*. The circles are the fit assuming 
a Gaussian P t (t) with the stress drop fixed to its real value. 
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Figure 9. Stress distribution estimated directly from the stress map shown in Figure 3b (circles), and inverted 
from the seismicity rate shown in Figure 8 (see model #4 in Table 1). The solid black line is the solution of equation 
(15)). The best fitting unconstrained Gaussian distribution is shown as a blue dashed line. The green dash-dot 
line shows the best fitting Gaussian with fixed stress drop. 
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Figure 10. (a) Seismicity rate for different values of the distance to the fault y/L decreasing from y/L = (top) to 
y/L = 0.2 (bottom), using the slip model shown in Figure 3a, and assuming Aa n = 1 MPa. (b) Standard deviation 
(solid line) and absolute value of the mean (dashed line, the average stress change is always negative) of the stress 
distribution as a function of the distance to the fault. 
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Figure 11. Variation of the effective Omori exponent with time, for a Gaussian stress distribution of mean 
—(To = — 3 MPa and standard deviation r* = 10 MPa. The dashed line is the exact solution (given by integrating 
numerically (6), and the solid line is the approximate analytical solution (26). 
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Figure 12. Seismicity rate for a Gaussian stress distribution (cr = —3 MPa, r* = 5 MPa), without (solid line) 
and with (crosses) heterogeneity of Aa n . Normal stress fluctuations are modeled by a log- normal distribution of 
average 1 MPa and standard deviation 7.3 MPa. The dashed line is a fit by the rate-and-state model (3), assuming 
Aa n is uniform. 



